//===========================================================================
/*!
 *  \file Potential.h
 */
//===========================================================================
#ifndef _Potential_H
#define _Potential_H


#include <math.h>

#include "OptCG.h"
#include "OptLBFGS.h"
#include "OptQNewton.h"
#include "OptFDNewton.h"
#include "NLF.h"
				 
#include "Individual.h"

using NEWMAT::ColumnVector;
using NEWMAT::Matrix;
using NEWMAT::SymmetricMatrix;

using namespace OPTPP;

enum POTENTIALMODEL{ OSS2, OSS2_ITER, TTM2F,TTM21F, TTM3F,
	Ghemical,    //Ghemical force field.
 MMFF94, //   MMFF94 force field.
 MMFF94s,   // MMFF94s force field.
 UFF,    // Universal Force Field
 TTM21F_ICE,
 TTM31F_ICE,
 QSPCFW_ICE,
  SPCF_ICE,
 HF_OSS2,
 EMPTY_EQU
};

				 
/*!
 *  \class Potential
 *  \brief Parent class of Potential models
 */
class Potential{	
	public:
		enum OPT_METHODS{ QNEWTON=0, CONJUGATEGRAD=1, LBFGS=2, FDNEWTON=3 , DFPMIN=4};
		Potential() { dim=0; model=0; unit=0; neval=0; 
			MaxStep=0.2; GradTol=1e-4; FcnTol=1e-6; nMaxEval=1000; opt_med=QNEWTON;
			debug_eval=0;
		}
		
		Potential(int _dim,int _model=0,int _unit=0) { SetDim_(_dim);	model=_model; unit=_unit; neval=0; 
			MaxStep=0.2; GradTol=1e-4; FcnTol=1e-6;	 nMaxEval=1000; opt_med=QNEWTON;
			debug_eval=0;
		}
		virtual ~Potential(){};
		
		
		virtual const char* classname(){ return "Potential";} //!< identifier
	

		virtual double evaluate_( double const *x ){};
		
		virtual void gradient_(double const* x, double* grad){
			if(!HasAnalyticalGradients()) gradient_num(x,grad);
		};
		
		double gradient_num(double const* xt, double* grad,double delta=1e-4);
				
		virtual void hessian_(double * _p, double**_gg,double delta=1e-4);
		
		double evaluate_(const ColumnVector& x ){	return evaluate_(x.data());	};
		
		void gradient_(const ColumnVector& x , ColumnVector& g ){
			gradient_(x.data(),g.data());
		};
		
		void gradient_num(const ColumnVector& x , ColumnVector& g ,double delta=1e-4){
			gradient_(x.data(),g.data());
		};
		
		virtual int optimize(double *x,int debug=0);
				
		virtual void InitInv(Individual& inv){};
		
		virtual double evaluate_(Individual& inv ){ inv.energy=evaluate_(inv.x);return inv.energy;};		
		
		virtual void gradient_(Individual&  inv){ gradient_(inv.x,inv.force); };		
		
		virtual double gradient_num(Individual& inv, double delta=1e-4){
			inv.rmsGrad=gradient_num(inv.x,inv.force);
			return inv.rmsGrad;
		}
		
		virtual int optimize(Individual& inv,int debug=0);
		
		//virtual void hessian_(Individual& inv, double**_gg,double delta=1e-4)
				
		virtual bool HasAnalyticalGradients() { return false; }
		virtual double ValidateGradients(double *x,int verbose=0);
		//virtual double eval_gradient_(double const* x, double* grad){};		

		int Model_(){ return model;} //!< get potential  model, see comm/def.h
		void SetModel_(int model_){ model=model_;}; //!< set potential model, see comm/def.h
		 
		int SetModel_(const char *e){ 			
			for(int i=0;i<EMPTY_EQU;i++){
				string t=Equation[i];
				if(t.compare(e)==0){ model=i; return model;}
			}
		};
		
		const char* GetUnit_(){
			switch(unit){
				case HARTREE: return "HARTREE";break;
				case KCALMOL: return "KCALMOL";break;
				case KJMOL: return "KJMOL";break;
				case DIMENSIONLESS: return "DIMENSIONLESS";break;
			}
		}  //!< return current unit of energy
		int Unit_() { return unit;}
		void SetUnit_(int unit_){ unit=unit_;}; //!< set unit, see comm/def.h
		
		int Dim_(){ return dim;};
		virtual void SetDim_(int dim_){ 		
			if(dim!=dim_){				
				init_x.Release();	final_x.Release();				
				dim=dim_;	init_x.ReSize(dim);	final_x.ReSize(dim);  
			}
		};
		
		virtual void SetExtra_(void *p){};
		
		virtual int ReadParam(const char *filename){}; //!<read parameters for potential model from file

		virtual int ReadParam(const double *alpha_){};
		
		int genome(double const* x);
		int genome(double const* x,int dim_);
		
		int MaxEval_(){ return nMaxEval;}
		void SetMaxEval_(int a){ nMaxEval=a;};
		
		double FcnTol_(){return FcnTol;}
		void SetFcnTol_(double a){ FcnTol=a;}
		
		double GradTol_(){return GradTol;}
		void SetGradTol_(double a){ GradTol=a;}
		
		double MaxStep_(){return MaxStep;}
		void SetMaxStep_(double a){ MaxStep=a;}
		
		int OptMethod_(){ return opt_med;}
		void SetOptMethod_(int a){ opt_med=a; }
		
		void SetDebugEval_(int debug=0){ debug_eval=debug;};
		
		int getNEvals(){ return neval;} 
		double getSeconds(){ return nsecs;} 
		void getFinalEval(double &fe){fe=final_eval;}
		void getFinalX(double *x){ for(int i=0;i<dim;i++) x[i]=final_x(i+1);}
		void getRmsGrad(double &fe){ fe=RmsGrad;}
		void getMaxGrad(double &fe){fe=MaxGrad;}		
		
		int checkConvg(){ return iConverged;}
		
		virtual void PrintInfo(){ cout<<"Base potetial"<<endl; };
		
		const char* equation(){ return Equation[model];};
		
		static const char* Equation[50];
		
		static double* AllocateArray(int x1){	double* d=new double[x1];	return d;}
		
		static double** AllocateArray(int x1,int x2){
			double** d=new double*[x1];
			for(int i=0;i<x1;i++) d[i]=new double[x2];	assert(d);
			return d;
		}
				
		static void FreeArray(double** d,int x1){
			if(d){	for(int i=0; i<x1; i++) if(d[i]) {	delete[] d[i]; d[i]=NULL;}
				delete[] d;			d=NULL;	
			}
		}
				
		static void FreeArray(double* d){	delete[] d; d=NULL;}
		
		static double ConvertUnit(double e,int from, int to){
			double ConvertTable[3][3]={1.0,HART_2_KCAL, HART_2_KJ,
				KCAL_2_HART, 1.0, KCAL_2_KJ,
				KJ_2_HART,KJ_2_KCAL,1.0};
			return e*ConvertTable[from][to];
		}

	protected:
		int opt_med; //!< optimization method
		int model; //!< index of model
		int unit; //!< index of unit system		
		int dim;	//!< number of dimensions 
			
		ColumnVector init_x;
		ColumnVector final_x;
		double final_eval;
		
		int nMaxEval;
		double FcnTol;
		double GradTol;
		double MaxStep;
				
		int iConverged;//!< check whether converged or not
		int neval; //!< number of evaluations performed
		double nsecs; //!< time(seconds) for optimization
		
		double MaxGrad,RmsGrad;
		
		int debug_eval;
	
	private:
		static Potential* pSelf;
		
		static void initgenome(int ndim, ColumnVector& x){				
			x=pSelf->init_x;
		}
		
		static void evalgenome_0(int ndim, const ColumnVector& x, double& fx, int& result){			
			fx  = pSelf->evaluate_(x);
			result = NLPFunction;
		}
		
		static void update_model(int, int, ColumnVector) {}
		
		static void evalgenome_1(int mode, int ndim, const ColumnVector& x, double& fx, ColumnVector& g, int& result){			
			if (mode & NLPFunction) {
				fx  = pSelf->evaluate_(x);
				result = NLPFunction;
			}
			if (mode & NLPGradient) {
				pSelf->gradient_(x,g);				
				result = NLPGradient;
			}			
		}		
		
		int optimize_QNewton(NLF1 &nlp);
		int optimize_CG(NLF1 &nlp);
		
		template<class T>	int optimize_nlf1(NLF1 &nlp);
		
		void lnsrch(int n, double xold[], double fold, double g[], double p[], double x[],double *fret, double stpmax, int *check);
		
		int dfpmin(double p[], int n, int iterMax, double *fret);
};

#endif